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An efRcient spectral interpolation routine for the TwoPunctures code 
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TwoPunctures is perhaps the most widely-adopted code for generating binary black hole "punc- 
ture" initial data and interpolating these (spectral) data onto evolution grids. In typical usage, 
the bulk of this code's run time is spent in its spectral interpolation routine. We announce a 
new publicly- available spectral interpolation routine that improves the performance of the original 
interpolation routine by a factor of <--T00, yielding results consistent with the original spectral inter- 
polation routine to roundoff precision. This note serves as a guide for installing this routine both in 
the original standalone TwoPunctures code and the Einstein Toolkit supported version of this code. 

PACS numbers: 04.25.dg, 04.40.Nr, 04.25.D-, 04.25.dk 



I. INTRODUCTION 



One of the pressing goals of numerical relativity is to calculate accurate gravitational waveforms from plausible 
astrophysical sources to help generate templates that will be used by ground based gravitational wave observatories 
such as LIGO [1,2], VIRGO [3, 4], TAMA [5, 6], GEO [7], KAGRA [8], and by proposed space-based interferometers 
such as cLISA/NGO [9] and DECIGO [10]. This task is far from trivial and the very first step in accomplishing it is 
CJ ', the generation of initial data that satisfy the Einstein constraints [11]. 

O" 1 To date there are multiple codes that solve the constraints of Einstein's theory of general relativity (see e.g. [12-17] 

5^ \ and references therein) . As the inspiral and merger of compact binaries such as binary black holes (BHBH) , binary 

bsO neutron stars and binary black hole-neutron stars, furnish some of the most promising astrophysical scenarios (both 

in terms of signal strength and event rates) for the generation of detectable gravitational waves, the solutions obtained 

^vq . using these initial data solvers are mainly focused on compact binary systems. 

£> " Most of the recent focus in gravitational wave template generation has been on binary black hole systems (see e.g. 

t**« . [18, 19]) and for this reason one of the most popular codes for generating BHBH initial data is the TwoPunctures code 

••/"") ■ [17], which has also been adopted by the publicly available Einstein Toolkit [20]. The popularity of this code stems 

from the fact that it is remarkably user-friendly, spectrally accurate and efficient in solving the Einstein constraints 

^— ' . ■ for BHBHs when both BHs are represented as punctures. Once the initial value problem has been solved, the initial 

data have to be mapped onto the dynamical evolution grids via interpolation. The TwoPunctures code offers two 

interpolation routines: i) a second-order polynomial interpolation routine, and ii) a spectral interpolation routine. 

The former is very fast but we have found empirically that it is not well-suited for dynamical evolutions. For this 

reason, it is the latter that is most widely used by us and most numerical relativity groups. 

In this brief note we provide a spectral interpolation routine for the TwoPunctures code (and installation instructions 
both for the Einstein Toolkit version of the code, i.e., the TwoPunctures thorn, and its standalone version) that is 
~ 100 times faster than the original spectral interpolation routine of the TwoPunctures code. This new routine saves 
many hours of computation, especially when high spectral resolutions are used. By no means do we claim that we 
have optimized the process, and one may find other ways of optimizing the performance of the code in general (see 
below). However, we find that the acceleration attained by our routine is sufficiently satisfactory, and we hope that 
by making this routine publicly available, the many numerical relativity groups that use TwoPunctures will benefit 
from this faster spectral interpolation routine. 

II. A FASTER SPECTRAL INTERPOLATION ROUTINE 

Our basic modification of the TwoPunctures code was stimulated by the observation that each time the original 
TwoPunctures code spectral interpolation routine is called it computes the spectral interpolation coefficients, given 
the values of the function at the collocation points, and then uses the spectral expansion to interpolate to any one 
Cartesian grid point. Typical high resolution evolution grids in finite difference codes employ 8-9 levels of refinement 
with resolutions of M/40 or higher, where M is the BHBH ADM mass. This amounts to a grid of about 10 5 — 10 6 zones 
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FIG. 1. Left Performance of the original TwoPunctures code interpolation routine normalized by the performance of our new 
spectral interpolation routine versus the total number of basis functions (or collocation points). Performance is measured 
in seconds of wall time. The inset shows the required wall time in seconds. The black curve corresponds to the original 
TwoPunctures code interpolation routine and the red curve to our spectral interpolation routine. The performance test is 
based on timing the interpolation operation for a grid consisting of 10 zones, after the code has solved for a particular BHBH 
configuration. The spectral resolutions used were (ni,n 2 ,n 3 ) = (10,20,20), (20,30,30), (20,40,40), (30,50,50), (40,60,60), 
(50,70,70). Right: Fractional difference between the interpolated values along the BHBH binary axis (the punctures lie at 
x/M = ±2.17) using the original TwoPuncture code interpolation routine and our new spectral interpolation routine. The 
results agree to machine precision. All runs were performed on a system with an Intel Core 2 Duo 6300 processor. The code 
was compiled with the Intel 11.1 CH — h compiler with -03 optimizations. Similar relative performance is found on other systems 
and with other compiler optimizations. 



for which interpolations must be performed. Calculating the spectral coefficients from the values at the collocation 
points is an expensive operation, which is exacerbated by computing them at every interpolated point. In our initial 
value calculations we use the TwoPunctures code with 50 x 70 2 basis functions, thus the original TwoPunctures code 
spectral interpolation routine calculates ~ 2.5 x 10 5 spectral coefficients each time it is called, i.e., ~ 10 6 times. This 
means that the spectral coefficients are recomputed ~ 10 11 times in this process. 

The total cost of this operation can be reduced significantly, if the spectral coefficients computation is a one-time 
operation. So, our new interpolation routine uses the Chebyshev and Fourier basis routines in the TwoPunctures code 
to compute the spectral coefficients once and for all, and then store them in an array. Another routine we developed 
takes the stored spectral coefficients as input and performs the spectral interpolation. To compute the required sums 
we use the partial summation method [21]. Our new routine still computes the values of the bases functions at the 
collocation points in order to perform the interpolation every time it is called, hence further optimization can take 
place by storing the values of the bases functions at the collocation points as well. However, we did not do so because 
the new interpolation routine performs sufficiently fast. 



In the left panel of Fig. 1 we show the performance of the original TwoPunctures spectral interpolation routine 
normalized by the performance of our routine against the total number of collocation points for a test case where 
we do 10 5 interpolations after solving for an initial BHBH configuration. It is clear that our routine accelerates the 
interpolation procedure. The higher the spectral resolution the greater the gain, and for resolutions of 50 x 70 2 we get 
more than a factor of 100 speed-up: our routine takes ~ 200s, while the original routine requires ~ 28000s to finish. 
We have also checked that the results with the two routines agree to machine precision as shown in the right panel of 
Fig. 1. 



III. INSTRUCTIONS FOR INSTALLING THE NEW INTERPOLATION ROUTINE 

Here we provide instructions for installing the new interpolation routine both in the TwoPunctures thorn of the 
Einstein Toolkit, and the standalone version of the TwoPunctures code. The routines for the Einstein Toolkit version 
of the TwoPunctures code can be downloaded from 

http://webusers.physics.illinois.edu/~vpaschal/TwoPuncturesET/. The routines for the standalone 
TwoPunctures code can be downloaded from 
http: //webusers .physics . illinois .edu/~vpaschal/TwoPunctures_Standalone/. 

A. TwoPunctures Einstein Toolkit 

We provide the routines 
void SpecCoef (int nl, int n2, int n3, int ivar, CCTK_REAL *v, CCTK_REAL *cf) 

CCTK_REAL PunctlntPolAtArbitPositionFast (int ivar, int nvar, int nl, 

int n2, int n3, derivs v, 
CCTK_REAL x, CCTK_REAL y, CCTK_REAL z) 

CCTK_REAL PunctEvalAtArbitPositionFast (CCTK_REAL *v, int ivar, CCTK_REAL A, 

CCTK_REAL B, CCTK_REAL phi, int nvar, 
int nl , int n2 , int n3) 

The first computes the spectral expansion coefficients (cf ) of a variable (v) given the values of variable v at 
the collocation points, where nl, n2, n3, are the number of basis functions in the A, B, (f> coordinates of the 
TwoPunctures code. The second routine applies spectral interpolation on the variable v to the desired point with 
Cartesian coordinates x, y, z, using the third routine which interpolates v using the corresponding TwoPunctures 
coordinates A, B,(f>. The following are the required steps to implement our new interpolation routine. 

1. For convenience we recommend that these routines be added in the file FuncAndJacobian.c of the TwoPunctures 
thorn. 

2. These routines must also be declared in the TwoPunctures. h header file. 

3. In the file TwoPuncture.c the following declaration statement must be added: 

static derivs cf_v; 

4. In the file TwoPuncture.c the following memory allocation call for storage of the spectral coefficients must be 
added: 

allocate_derivs (&cf_v, ntotal) ; 

5. Following memory allocation, the new variables should be initialized (e.g.): 

for (int j = 0; j < ntotal; j++) 
{ 

cf_v.dO[j] = 0.0; cf_v.dl[j] = 0.0; 

cf_v.d2[j] =0.0; cf_v.d3[j] =0.0; 

cf_v.dll[j] = 0.0; cf_v.dl2[j] = 0.0 

cf_v.dl3[j] =0.0; cf_v.d22[j] =0.0 

cf_v.d23[j] =0.0; cf_v.d33[j] =0.0 
} 

Note that there is a redundancy, as cf _v . dO is the only required variable, but as the memory footprint of the 
cf _v derivs struct is small this redundancy is of minor significance. 

6. In the file TwoPuncture.c following the function call 



F_of_v (cctkGH, nvar, nl, n2, n3, v, F, u) ; 

which follows the loop checking for convergence of the puncture masses, the call to the calculation of the spectral 
coefficients should be made 

SpecCoef (nl, n2, n3, 0, v.dO, cf_v.d0); 

7. All calls to PunctlntPolAtArbitPosition should be replaced by corresponding calls to 
PunctlntPolAtArbitPositionFast, but the most important change is to replace the call 

U = PunctIntPolAtArbitPosition(0, nvar, nl, n2, n3, v, xl, yl, zl) ; 

by the call 

U = PunctlntPolAtArbitPositionFast (0, nvar, nl, n2, n3, cf_v, xl, yl, zl) ; 

8. Memory for the spectral coefficients should be freed at the end of the TwoPuncture.c file 
free_derivs (&cf_v, ntotal) ; 

B. Standalone TwoPunctures code 
For the standalone TwoPunctures code we provide the following routines: 

void SpecCoef (parameters par, int ivar, double *v, double *cf) 

double Spec_IntPolFast (parameters par, int ivar, double *v, double x, double y, double z) 

double Spec_IntPolABphiFast (parameters par, double *v, int ivar, double A, double B, double phi) 

These routines do precisely the same calculations as the routines SpecCoef , PunctlntPolAtArbitPositionFast , 
PunctEvalAtArbitPositionFast in the TwoPuncturesET version of the code. 

The steps to implement our new interpolation routine in the standalone version are similar to the ones we outlined 
in the previous section III A 

1. These routines can be added to the file FuncAndJacobian.C of the TwoPunctures code. 

2. These routines must be declared in the TwoPunctures. h header file. 

3. In the file TwoPunctures. C the following declaration statement must be added: 

derivs cf_v; 

4. In the file TwoPunctures. C memory should be allocated for storage of the spectral coefficients: 
allocate_derivs (&cf_v, ntotal); 

5. In the file TwoPunctures. C following the call 
Newton (par , v) ; 

the calculation of the spectral coefficients should be performed: 
SpecCoef (nl, n2, n3, 0, v.dO, cf_v.d0); 

6. Calls to 

Spec_IntPol(par , 0, v.dO, x, y, z) ; 

should be replaced by calls to 

Spec_IntPolFast (par , 0, cf_v.d0, x, y, z) ; 

7. Memory for the spectral coefficients should be freed at the end of the TwoPunctures. C file 
free_derivs (&cf_v, ntotal); 
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